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I. INTRODUCTION 



The effects of mountain ranges on synoptic scale motions 
have been shown (Reiter, 1971) to be quite influential. Yet, 
the physical cause and mathematical description of these ef- 
fects are poorly understood. The large-scale flow over and 
around mountain ranges is not properly observed or mathema- 
tically well modeled. This was quite evident to this author, 
having spent three years of forecasting weather for the Euro- 
pean-Mediterranean area, where the Alps, Pyrenees, and other 
mountain ranges have a marked meteorological effect. Due to 
the grid size of most primitive equation models, a nested, 
finer resolution grid is required to depict accurately the 
flow over and around mountains. This is attempted in this 
study . 

The general circulation model used in this study is a 
modified version (Monaco, 1975) of the UCLA global predic- 
tion model. To reduce computer memory and run-time require- 
ments prior to nesting a finer resolution grid, a channeling 
of the global model was undertaken. Procedures had to be 
developed to take into account the scheme C (Monaco, 1975) 
grid arrangement of the model and the higher order form of 
the flux terms in the primitive equations. 

In nesting a fine resolution scheme C grid, certain con- 
ditions had to be considered to take into account the stag- 
gered grid arrangement. Unlike the majority of nested grid 



11 



models presently available which use unstaggered grid ar- 
rangements, the indexing scheme used for each variable must 
be considered in determining interface boundary conditions. 

The solutions obtained are divided into three categories: 
a) coarse mesh solution, b) uniform fine mesh solution, and 
c) nested solution. Comparisons are made between solutions 
with no mountain present, and then with a mountain. It is 
shown that a number of problems remain for the type of 
channeling and nesting used in the present model. 
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II. 



MODEL DESCRIPTION AND 
INITIALIZATION PROCEDURES 



The general circulation model used in this research has 
been described by Monaco and Williams (1975). It is a modi- 
fied version of the UCLA global prediction model detailed by 
Arakawa and Mintz (1974) . Our model assumes the atmosphere 
is adiabatic and frictionless and contains no sink or source 
terms. Numerical integration is carried out on a staggered, 
spherical, sigma coordinate system. 

A. VERTICAL COORDINATE 

The vertical coordinate used in the model is the non- 
dimensional a-coordinate, fig. 1. It is defined as follows: 




where tt is the terrain pressure. The pressure is given by 
p , the surface pressure by p , and the constant tropopause 
pressure by p . The vertical boundaries of the coordinate 
system, the earth's surface and the constant tropopause 
pressure level, follow from (2.1), such that: 

a = 0 at p = p t (= 200 mb) 

a = 1. at p = p s . 
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The boundary condition at the tropopause pressure level, 

H rr 

and at the earth's surface, is: 6 (= ^r) = 0 . 

B. HORIZONTAL DISTRIBUTION OF VARIABLES 

The dispersion of local accumulations of wave energy 
generated by gravity-inertia waves is an important process 
in any atmospheric numerical model. Winninghoff (1968) 
showed how this geostrophic adjustment depended on the dis- 
tribution of variables over the grid points. Five possible 
distributions of the dependent variables <p , T, v, and u are 
shown in fig. 2. The original UCLA general circulation 
model, described by Arakawa and Mintz (1974), used scheme B 
for the horizontal distribution of variables; this research 
model uses scheme C. 

Winninghoff showed that schemes B and C, in a one-dimen- 
sional case, adequately simulated geostrophic adjustment. 
However, in a two dimensional case, scheme C was shown to 
be the best lattice to simulate the geostrophic adjustment 
process. In using scheme C, though, two conditions can cause 
problems. The first is due to the two velocity components 
not being carried at the same latitude, causing some diffi- 
culty with Coriolis force terms. The assumption of a con- 
stant Coriolis parameter could eliminate this problem. The 
second problem occurs if A^y^ is less than or close to one. 
The quantity A R , the Rossby radius of deformation, is given 
by 



A r = /p /£ 
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where H is Che mean value of h , Che depch of Che acmos- 
phere. The siCuacion where £ 1 is rare and, chere- 

fore, scheme C is used. 

C. PRIMITIVE EQUATIONS 

The primicive equacions are wriCCen using Che orchogonal 
curvilinear coordinaCes , £ and n , where £ = X (longi- 

Cude) and n = d (laCiCude) . Scaling faccors, m and n , 
are defined as follows: 



1 x J 1 

— = a cos d> and — 
m n 



a . 



The horizoncal velociCy componenCs, Cemperacure, and surface 
pressure, Che prognoscic variables, are governed by Che fol- 
lowing primicive equacions: 



3 / TT \ , 3 /ITU \ , 3 /1TV \ 

scW u) + -5I ( — u) + ^ ( -5T u) 



+ 



£_(££ u ) 
3a J 



( 2 . 2 ) 



on 



j. / 3 1 

+(v 3£ n 



3 L , , tt 9 4) 

u aTT a )|wv + H [ 3f 



3 TT 

aa tt?-] 



35 ' 



JL 

mn 



V ) + v) + v ) + l-fli v) 

3 C mn v; + 3£^n v; 3n^m v ' 3aW v; 

, r f , / 3 1 ,, 3 1 n 1w „ , TT 9<j) , 3 tt _ 

mn 9^ n 9n nr m 9n 9n 



( 2 . 3 ) 



— F 
mn t) ’ 



3 / 77 m\ i 3 / TT U m \ i 3 / TT V rp \ 

C P T) + T5 ( tt c P l)+ 37T ( -ir C P T) 



+ (-E-) K c e) 

^p ' 3a W p ; 



( 2 . 4 ) 



.3 / tt \ , u 3 tt,v 3tt , n n 

" aa ‘itW + 5jT + i^' + 75 Q 
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and 




(2.5) 




( 2 . 6 ) 



6$ = -ira6a . 



(2.7) 



A complete listing of the symbols used in the preceding equa- 
tions may be found in the front of this report. 

The vertical a-velocity is given by: 



In determining geopotentials, conservation of total energy 
under adiabatic and frictionless processes, and conservation 
of 9 and ln0 , integrated over the entire mass under adia- 
batic processes, are needed. In addition, the vertical dif- 
ference scheme must maintain the property of the vertically 
integrated horizontal pressure gradient force. This is a- 
chieved essentially by combining an alternate form of the 
hydrostatic equation, written in terms of potential tempera- 
ture (2.9), with the vertical integral of hydrostatic 
equation (2. 10) : 



o 




( 2 . 8 ) 



o 



3$ = c 




(2.9) 




$) da . 



( 2 . 10 ) 



o 



o 
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The system is complete when the above two diagnostic equa- 
tions are included with the primitive equations. 



D. ANALYTIC INITIALIZATION 

The initial wind and pressure (geopotential) fields are 
defined mathematically in this model which has some advan- 
tages over obtaining the fields from a weather map. The 
analytic initial conditions: 1) simplify the task of inter- 

polating and balancing initial fields from constant pressure 
surfaces to sigma surfaces, 2) minimize imbalances between 
wind and pressure fields thereby reducing gravity-inertia 
waves . 

Kaurwitz (1940), in defining the initial velocity fields, 
solved the linearized, non-divergent vorticity equation for 
the stream function. The initial distribution of geopoten- 
tial was obtained by using this solution in the forcing func- 
tion of the non-linear balance equation (Phillips, 1959). 

Initialization of the model, based on the aforementioned 
solutions, was achieved using the following equations: 



2 

= -Ba sm <}) 



= a 2 A (<f>) 



u 



= _ \ = ® a cos $ 

a 9(j) 



( 2 . 11 ) 



v 



a cos <j> 3 A 



1 3 4 * __ q 

n o A 1 1 U • 



where 
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These equations represent solid rotation of the atmosphere 
because of the absence of vertical wind shear. The horizon- 
tal or zonal wind distribution varies as the cosine of the 
latitude and the magnitude of the constant, B. The zonal 
wind distribution in the model varied from u = 15.2 m-sec’^ 
at 34°N to u = 13.8 m-sec ^ at 58°N. The initial tempera- 
ture field was determined in accordance to the NACA standard 
atmosphere as defined by Haltiner and Martin (1957) . 



E. PROCEDURE FOR INTRODUCING TOPOGRAPHY INTO MODEL 

The scheme used to introduce mountain topography into 
the model is described in detail by Hayes (1977). Initially, 
heights at the surface were set equal to zero everywhere. 
During the model integration, the heights of the mountain 
grid points were increased as a function of time until the 
desired heights were achieved. The problem with wind field 
initialization in the vicinity of the mountain precluded 
starting with a mountain "already built". The following 
equations were used at the mountain grid points: 



,# • 2 / TT t \ . 

= M sin (jt~' * t - t p 



sf c 



( 2 . 12 ) 



= M 



t > t 



where $ s f c is t ^ ie geopotential height at time t , M is 
the final geopotential mountain height, and t is the 
period of mountain "building". Values used in this research 
are : 

^sfc(max) * 7500 • 

t =12 hours . 

P 
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III. CHANNELING THE COARSE GRID MODEL 



In any numerical modeling, considerations of computer 
run-times and memory requirements are necessary. This is 
especially true when developing and testing fine mesh grids 
where the grid spacing, and therefore the time step, are 
small. With this in mind, an attempt to channel the global 
prediction model (Arakawa, 1972; Arakawa and Mintz, 1974; 
Monaco and Williams, 1975) in the north-south direction was 
made prior to nesting of a finer resolution grid. 

A. PROCEDURES FOR CHANNELING 

In the initial phase of this research, the assumption 
of no variation in the north-south (n) direction was con- 
sidered. This led to reducing the amount of north-south 
pressure grid points in the model from 46 to 7. In accom- 
plishing this, the pole modification schemes for the 
various predictive equations were eliminated. 

In addition to removing the smoothing effect of the pole 
modifications, the effect of another scheme which zonally 
smooths the £ -components of the pressure gradient force and 
the divergence was reduced. This zonal smoothing eliminates 
computational instabilities due to the convergence of 
meridians at the pole. This smoothing is a function of lati- 
tude, grid size, and wave number; therefore, due to the 
channel being centered at 46°N with a width of only 8° lati- 
tude, its effect was much reduced. 
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Although the above routines might have had a damping 
influence on some of the large variations encountered with 
the model, the elimination or reduction of these routines 
did not generate any oscillations. Rather, these oscilla- 
tions were generated by the channel boundary conditions. 

B. BOUNDARY CONDITIONS 

Simple cyclic conditions east-west were used, but the 
north- south boundary conditions proved to be quite formid- 
able. The main constraint, for the north-south direction, 
was that of no net mass flux in the q-direction. The two 
main problems in this regard were the balancing of pressure 
gradient term to the u-velocity component and the integrated 
mass flux in the q -direction. Compounding the problem fur- 
ther were the scheme C grid arrangement and the complex flux 
forms of the advective terms in the equation of motion. 

The scheme C grid point did not allow easy selection of 
the boundary walls, due to the asymmetry of the variables 
at each (I,J) grid point. The flux terms of the advective 
equations involved a second order Jacobian scheme which 
required specification of some terms outside the walls. 

These problems were implicitly inherent in all boundary con- 
siderations . 

For the constraint of no mass flux in the north- south 
direction to be achieved, the v-equation of motion, (2-3) , 
integrated in the vertical and ^-directions must balance. 

The three terms that must be specifically dealt with along 
the boundaries are the integrated pressure gradient, normal 
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component of mass flux, and Coriolis terms. (Details are 
given in Appendix A.) 

The method used to prevent net mass flux across the 
walls was to reflect the normal component of the mass flux 
from the immediate interior grid point by one opposite in 
sign but equal in magnitude just outside the wall. This was 
accomplished by the following equations: 



(7rv ¥>i,2 + (itv ¥>i,i = 0 



m ^ I , Jmax + (Trv ra^I.Jmax+l 0 



(3-1) 

(3-2) 



The term (ttv ^-) ^ j nax+ ]_ is a mass flux quantity outside 
the boundary wall, which was needed in the advective terns 
of the equations of motion. 

Although this procedure assured no mass flux across each 
boundary, it did not assure zero net mass flux over the 
entire field. This is due to the inequality between the 
northern and southern mass flux terms, i.e.: 

TfA.l + (uv TTh.Jmax+l + 0 • (3 ' 3) 



This imbalance did contribute to total field mass variations 
with time. 

Along the walls the balance between the pressure gradi- 
ent and the geos trophic u-component, 



f TT 
mn 



u 



tt r 3b , 3ir , 
+ aa 3TT ] 



(3-4) 
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was difficult to maintain. The gradient of terrain pres- 

3 7T 

sure, , was restored to its initial value after each 

time step. This allowed the pressure to. vary along the walls, 
but maintained the initial pressure gradient. The geopoten- 
tial height gradient was controlled in a different way. 

Since geopotential height fields are not stored but calculated 
at every time step, restoring the gradient -1^ to its initial 



3n 



value was not straightf orward . 



3 $ 



The method used to specify was through the parameters 

that $ was a function of, i.e. temperature, pressure, and, 
indirectly, the velocity components. Since the boundary grid 
points were not predictive points, the values needed at the 
boundary walls were provided by extrapolation schemes using 
interior predictive values. Since = 0 at a = 1 , but 



3tt 



increases with elevation as act decreases, any imbalance 

in the parameters used to control 9 become more significant 
with height. With or without a mountain in the model, the 
largest variation of the v-component always occurred at sigma 
level 1 and decreased downward. This was evidently due to 

r\ i 

the large term in the v-component equation. 

With the inequality of normal mass flux terms plus the 
imbalance in the boundary zonal flow, there was a variation 
in the total field mass, and the oscillations were evident 
in the surface pressure fields. This variation of pressure 
along the north-south direction with time is shown in Fig. 3. 
The approximate period of the oscillations is 20 hours, which 
indicates that the forcing might originate from the southern 
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boundary where the inertial period is 21.5 hours. This 
might indicate a net mass inflow from the south in the early 
stages . 

One attempt was made to eliminate these pressure, i.e. 
mass, variations. The change in terrain pressure between 
two time steps, averaged over the entire grid, was subtracted 
from the forecast pressure fields. As shown in Fig. 3, this 
reduced the oscillations substantially but, as shown in 
Fig. 4, caused an increase in the temperature variation. This 
temperature increase then caused problems in the calculation 
of geopotential heights along the boundaries. The procedure 
was discarded, thereby allowing net mass changes over the 
entire field to occur. 

The phase and period of these pressure oscillations were 
dependent on the grid and domain size. The oscillations in 
a model with a grid distance of 4° latitude versus 1.5° lati- 
tude are compared in Fig. 5. In addition, the domain size 
of the 1.5° latitude grid is one- third that of the 4° latitude 
grid. The phase shift of six hours between the two grid 
spacing is obvious, but there is also a slight period change. 
The smaller grid has a period of approximately 18 hours. 

This is the inertial period for the southernmost run which 
adds evidence that the mass variations are forced from the 
south. 

C. GRID ARRANGEMENT AND ADVECTIVE TERMS 

As mentioned in the preceding section, the scheme C grid 
arrangement and the complex form of the advective terms 
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compounded Che problems at the boundary walls. The follow- 
ing is a brief description of these problems. 

1. Scheme C Grid Arrangement 

In Figs. 6A and 6B the parameters enclosed by a tri- 
angle comprise a scheme C (I,J) grid point. This indexing 
nomenclature lends itself readily to placing a boundary wall 
not at one latitude but two. This was used in this research, 
and proved to cause more problems than it solved. 

Since the boundary walls are not predictive points, 
values of tt , T , U , V must be provided to allow time 
stepping of the interior points. As discussed in section B, 
the extrapolated methods used in this research were not com- 
pletely satisfactory. An example of this is the divergence 

calculation at ir-point (1+1, Jmax-1) . The j— - ’’- -- term is 

o t, 

based completely on predictive values, where the ^ ^ TTV -- 

term is not. Since the same gradient of tt and V is main- 

tained between (Jmax-1, Jmax-2) and (Jmax, Jmax-1), an 
artificiality could develop with the presence of a distur- 
bance along the boundary. As long as an indexed (I,J) grid 
point is used as a channel wall, this problem will exist. 

An alternative channel wall boundary was considered, 
but not tried. It was to have had an unequal number of 
rows between the parameters U , <J> , T , ir , d and V in 

the north-south direction. This would allow the U , <p , 

T , tt latitude row to be predictive along both the north 
and the south boundaries. This would also eliminate the 
problems of imbalances caused by the use of extrapolation 
methods for providing values at the walls. Indexing does 
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become much more complicated, though, and this was the pri- 
mary reason it was not attempted. 

2 . Flux Form of the Advective Terms 

The prevention of non-linear instabilities in a non- 
divergent flow can be achieved by conserving the mean square 
vorticity, and mean kinetic energy with time. The global 
model used in this thesis was written to simulate slowly 
changing quasi-geostrophic motion. Prevention of non-linear 
instability called for a second order finite difference 
Jacobian. This second-order Jacobian was used by Arakawa 
(1974) to write the non-linear advective terms of the primi- 
tive equations. Due to its second-order nature, the mass 
flux expression for a specific grid point required advective 
mass flux terms as far away as two (I,J) grid points. This 
presents a problem as the boundary is approached by requiring 
the specification of extra mass flux terms outside of the 
channel wall (line Jmax+1 in Fig. 6A) . 

The value of the normal mass flux outside the channel 
wall is a reflected value of Jraax. The mass flux at Jmax, 
though, is an extrapolated value and not a predictive value. 
This is due to the scheme C grid arrangement and choice of 
location of the channel wall. At the southern channel wall, 
a predictive mass flux quantity at J=2 is reflected out to 
J=l. 

If the channel wall is placed at the latitude row 
having <p , T , it , U , and b , then a simplified flux scheme 
would be required. Although no extra flux terms would be 
needed for predicting temperature, another external (rru 
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term would be needed at J=0 and J=Jmax+l to calculate momen- 
tum flux centered on a U point (Fig. 6B) . It would probably 
be adequate to simply set the following: 
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To briefly summarize, the scheme C grid arrangement 
and indexing readily leads to certain boundary wall condi- 
tions. These proved to be unsatisfactory, due to extrapola- 
tion techniques for providing the channel wall values. 
Another channel wall scheme was presented which would 
eliminate extrapolating values to the wall. However, with 
the integrated mass flux in the ri direction not being zero 
with the changing of the flux scheme for the advective terms 
and with the large indexing problems, this scheme can still 
create difficulties. Time did not permit attempting the 



latter scheme. 
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IV. NESTED MODEL 



A. SCHEME C NESTING 

The nested grid used in this thesis consisted of two 
uniform grids with the fine mesh grid centered in the chan- 
neled coarse mesh grid (Fig. 7). Due to the scheme C stag- 
gered grid arrangement, the ratio of the coarse mesh length 
to the fine mesh length is three to one, which corresponds 
to a reduction of latitude spacing of 444 km to 148 km at 
the center latitude. The longitude reduction is from 348 km 
to 116 km, also at the center latitude. The total number of 
grid points in the coarse mesh for one sigma level is 20 x 7, 
while in the fine mesh it is 19 x 7. 

Although the fine mesh grid is placed at the center of 
the coarse mesh, there is a minimum number of J rows required 
to the north and south of the nested grid to allow inter- 
action between the grids at predictive points. Due to the 
extrapolation schemes used at the coarse mesh boundaries, at 
least two J rows are needed to the north and south of the 
nested grid. If the channel walls are changed to consist 
only of the <j> , tt , T , U latitude row, then this require- 
ment would be unnecessary. 

The physics of the model for the nested fine mesh are 
identical to those of the coarse mesh domain. The only 
difference between grids was the complete elimination of the 
zonal smoothing routine used in the coarse mesh. 
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B. BOUNDARY CONDITIONS 

There are numerous approaches for specifying the inter- 
face boundary conditions in a nested grid model (Ookochi, 
1972; Harrison and Elsberry, 1972; Chen, 1973; Moss and 
Jones, 1973; Chen and Miyakoda, 1974; Jones, 1976; Madala, 
1977; and Miyakoda and Rosati, 1977). They range from in- 
terpolation schemes to using the finite difference equations 
in providing the various boundary values. The method used 
in this paper was the Dirichlet condition with local bound- 
ary smoothing (Chen and Miyakoda, 1974). 

The boundary points are provided all necessary values 
of U , V , 7T , and T for both one and two way interaction 
approaches from the solutions of the coarse mesh grid. 

These boundary values vary with time as the coarse mesh ad- 
vances in time. Linear interpolations in time and space are 
provided for the fine grid where necessary. Due to the flux 
form of the advection terms, additional mass flux terms 
(ttv and (ttu — ) are required outside the interface. 

These values are also interpolated from the coarse mesh. 

As described by Chen (1973) and Chen and Miyakoda (1974) , 
this method of providing all variables at the boundary is an 
over specif icat ion . This can cause a computational mode to 
develop at the boundaries. A method to reduce this computa- 
tional mode, but not to suppress significantly the physical 
mode is to apply a local smoother to all variables at the 
grid points next to the boundaries. Chen and Miyakoda (1974) 
suggested the use of a simple filter (1-2-1) applied normal 
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to the boundary at the adjacent row. In this research a five 
point filter was used to help suppress lateral waves along 
the boundaries. For example, the smoothing routine for the 
U-component at the west boundary is as follows: 



U 2, j ( smooth ed) 7 u 2,j + F^ U 1, J +U 2, J +U 2, J+l +U 2, J-l^ * 



This filter was not applied after every time step, but after 
every 6-3-3-3 time steps (Fig. 3). 

As with the coarse mesh boundary conditions , the scheme 
C grid arrangement and flux forms of the advective terms af- 
fected the choice of the interface walls. A scheme C (I,J) 
grid point was used for the interface walls. In Fig. 7 the 
points surrounded by the box are predictive points, whereas 
fine mesh lines, I, J, 1+9, J+6 are interface, non-predictive 
points . 

Unlike the channeled coarse mesh grid, there are no con- 
servation of mass or energy requirements for the nested fine 
mesh. Therefore, free exchange across the interface is 
necessary. Providing the boundary values with local smooth- 
ing was insufficient to maintain free exchange of mass and 
energy across the interfaces. This can be seen in Fig. 15b, 
which shows the reflections of the U-component from the east 
interface. A method suggested by Madala (1977), which con- 
sisted of calculating the mass convergence of the coarse mesh 
at the interface to provide values for certain boundary fine 
mesh mass flux terms, was attempted without success. (Refer 
to Appendix B for details of the divergence boundary condition.) 
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Problems arose at the west and north interface walls where 
the horizontal flux terms were removed considerably from 
the interface. Changing the interface walls to consist only 
of the tt , T , and <j> row would have allowed better conser- 
vation of divergence between grids. This change, though, 
would have entailed changing the flux forms of the advective 
terms and considerable indexing changes. Time did not permit 
complete testing of this boundary condition. 

C. MARCHING PROCESS FOR THE NESTED MODEL 

A schematic of the Matsuno-leapfrog time integration 
scheme for both the coarse and fine mesh grids is shown in 
Fig. 3. Each sequence is 45 minutes in length. This time 
sequence applies the Matsuno scheme less often in the fine 
mesh, which results in less damping of the high frequency 
oscillations. This is in contrast to Chen and Miyakoda 
(1974), who used the Euler-backward time step solely, and 
Miyakoda and Rosati (1977) , who used a time filter (Asselin, 
1972) after each time step. 

Due to the oscillations generated at the boundaries of 
the coarse and fine meshes, it was difficult to evaluate the 
above time stepping scheme. Computational modes did not 
appear to be significant in the various model cases studied. 
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V. 



DESCRIPTION OF THE 
EXPERIMENT AND RESULTS 



The main objective in all nested grid models is to show 
that the nested solution gives results similar to those of 
a uniform fine mesh. This applies to situations in which 
the resolution of the coarse mesh grid is inadequate and the 
fine mesh resolution is sufficient. An attempt was made in 
this research to obtain a nested fine mesh solution for flow 
over mountains. 

The following is a listing of the three basic models 
studied in this research. 

(1) Coarse mesh solution: 

a) without a mountain 

b) with mountain A 

c) with mountain B 

(2) Uniform fine mesh solution: 

a) without a mountain 

b) with mountain A 

(3) Nested solution: 

a) without a mountain, one way 

b) without a mountain, two way 

c) with mountain A, one way 

d) with mountain A, two way 

Mountain "A" is the case in which the mountain is built com- 
pletely across the channel, whereas mountain "B" is built 
only half-way across the channel. Results to 12 hours are 
compared, due to problems developing after that time. 
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A. COARSE MESH SOLUTION 

Hayes (1977) studied the flow over mountains using the 
sane primitive equation model (Monaco, 1975), except he 
used the full global model and not a channeled version. Some 
of his results which can be used in comparison are as follows: 
(a) the mountains interact with the atmospheric flow to in- 
duce ridging over the mountains and troughing on the down- 
stream side; (b) damping of small scale features with height; 
and (c) increasing the grid resolution causes the scale of 
the velocity component response to decrease, but the magnitude 
remains about the same when compared to the coarse grid 
resolution. The grid size that Hayes used for the coarse 
mesh, 4.5° long by 4.0° latitude, is the same as used by 
this author. 

Prior to testing the channeled model with a mountain, a 
case was computed without a mountain. This was done to check 
if the initialization scheme, which used the linear balance 
equation, was satisfactory for the channeled model. The 
model was run to 56 hours with the maximum variation in both 
velocity components not exceeding + 1 m/sec over the initial 
values. Some adjustment was expected because the finite 
difference form of the analytic solution is inexact. The 
terrain pressure does experience a variation with time (Fig. 
3). This is due, as explained in Chapter III, to the inte- 
grated mass flux in the north-south direction not being equal 
to zero. Therefore, other than the pressure variation, the 
model does appear to be in good balance. 
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Next, a meridional mountain range 18° x 750 m (width by 
height) extending across the entire channel was considered. 
Since the distance between grid points in the ^-direction is 
4.5°, this allows only three grid points at each latitude to 
represent the mountain. The middle channel U and V compo- 
nents after six hours are shown in Figs. 9A and 9B for three 
levels. The mountain has developed to half of its height by 
this time. Upstream, the flow is deflected northward by the 
mountain and downstream the flow is deflected to the south. 

At some distance from the mountain both upstream and down- 
stream, variations due to the mountain are less evident. 

Figure 9B shows an increase in the U-component in the vicinity 
of the mountain top. The increase can partly be explained by 
the use of the conservation of potential vorticity equation, 

C + f 

— = constant , 

where z, is relative vorticity, f is the Coriolis para- 
meter and Ap is the pressure depth of a column. In order 
to conserve potential vorticity on the upstream side the flow 
turns anticyclonic and therefore an increase in the V-compo- 
nent results. Over the mountain crest the flow turns cyclonic 
and the V-component decreases, but due to the decreased column 
depth the conservation of potential vorticity is achieved by 
an increase in the U-component. 

The vertical velocity field, shown in Fig. 10, indicates 
rising motion on the upstream side with sinking motion just 
over the mountain crest and downslope. In addition, a 
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secondary region of rising motion is indicated downstream. 
This can be indicative of a mountain lee trough. 

The temperature patterns for the middle of the channel 
at the six different levels is shown in Fig. 11. Note the 
cooling on the ups lope side of the mountain and the warming 
downstream. This is a direct result of the vertical veloc- 
ities induced by the mountain. It must be noted, though, 
that the temperature field over the entire channel increased 
on the average of 10 °K at all levels from initial tempera- 
tures. This was due in part to term (1) of the equation in 
Appendix A, the integrated n-component of momentum and the 
imbalance of mass flux terms between the north and south 
boundaries . 

At 12 hours, the mountain has attained its final height. 
The horizontal wind profiles for this time are shown in 
Figs. 12A and B. The U-profiles appear reasonable with the 
highest zonal velocity over the mountain crest, and also a 
decreasing magnitude with height. The V-profiles, though, 
begin to show problems. At level six the flow is still being 
deflected in the proper way, but problems are appearing at 
levels one and three. Although they, too, are showing 
proper deflection, the maxima of level one is the same magni- 
tude as level six. This is unexpected, based on Hayes' 
results showing dampening with height. 

The large unexpected maxima at level one is possibly due 
to two factors. The first is due to the method of specifying 
3 — at the boundaries (refer to Chapter III for details.) 
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The second is the apparent development of a meridional, ver- 
tical oscillation as a result of the vertical velocity pat- 
tern induced by the mountains. Note in Fig. 12A the apparent 
compensating northward flow downstream at level six in re- 
sponse to the southward flow at levels one and three. In 
Fig. 13 the largest vertical velocities at level five, and, 
therefore, vertical transport of momentum, are shown to be 
occurring along the channel walls in the vicinity of the 
mountain. This pattern of the largest vertical velocities 
occurring at the channel wall-mountain intersection con- 
tinues with time and dominates the vertical velocity patterns 
after 24 hours. 

To estimate the effect of the mountain-channel wall in- 
tersection, a test case was computed using the half -channel 
mountain "B". The horizontal velocity profiles were similar 
to the channeled mountain, except for a slight increase in 
the zonal flow north of the mountain. The 12 hour vertical 
velocity pattern, Fig. 14, associated with level five of 
mountain "B" is slightly different, especially along the 
north wall, when compared to the pattern of mountain "A", 

Fig. 13. At 12 hours the largest vertical velocities appeared 
along both the north and the south walls. Even with the 
mountain not intersecting the north wall, vertical velocities 
of the same magnitude as in Fig. 13 developed. Further tests 
would be required to determine the minimum distance needed 
to prevent interaction of the flow around the mountain and 
the channel wall. 
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B. UNIFORM FINE MESH 



The sane geographic domain as in the coarse mesh model 
was used in the fine mesh model with one- third the grid dis- 
tance. Due to the increase in the computer memory require- 
ments (9X) and CPU (Central Processing Unit) time (27X) , 
associated with this solution, few case studies were com- 
pleted. 

The case in which there was no mountain was computed to 
evaluate the initial balance and pressure, i.e. mass, varia- 
tion over the total field. Results through 12 hours showed 
maximum velocity variations of + 1 m/sec. The average pres- 
sure over the domain did vary with increases of . 3 mb at six 
hours and 2.6 mb at 12 hours. The phase and period of these 
pressure oscillations were difficult to calculate as the re- 
quired long computer runs had to be limited due to the exces- 
sive processing times. It appears, though, to have a phase 
close to that of the coarse grids phase, but a smaller 
magnitude . 

Regrettably, due to time limitations, the uniform fine 
mesh with a channel mountain was not run for a sufficiently 
long time to use in this study. Results achieved by Hayes 
(1977) in increasing the horizontal grid resolution (£- 
direction only) can be used as a basis of what to expect. 
Mainly, he found that by increasing the resolution in the 
^-direction the velocity component response decreased, but 
the magnitude of that response remained approximately the 
same when compared to the coarser resolution grid. 
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C. NESTED SOLUTION 

Two methods of interaction are employed with the nested 
grid model. The first method is one-way interaction; values 
from the coarse mesh are supplied to the fine mesh boundary 
at appropriate times, with no corresponding changes of coarse 
mesh interior values from the nested fine mesh. The other 
method is two-way interaction. Although boundary values for 
the fine mesh are still provided by the coarse mesh, replace- 
ment of coarse mesh values by interior nested fine mesh 
values at coincident points is done at appropriate times. 

For example: v(I+3,J+2) of fine mesh replaces v(I+2, J+2) 

of coarse mesh. A five point averaging scheme surrounding 
each fine mesh coincident point is used to calculate the 
replacement fine mesh value. 

Prior to testing the nested grid with a mountain present, 
one case with each type of interaction was run without a 
mountain to check the initial balance. As mentioned in sec- 
tion 1 of this chapter, the coarse mesh is in good balance. 
The v and u components in the middle region of the fine mesh 
after six hours are shown in Figs. 15A and B. The corre- 
sponding coarse mesh solution can be seen at either end of 
each plot. The v-component profile shows level l(a = .03) 
to have the largest value. This can be due to the following: 

1) the oscillations of the height surfaces in the north- 
south direction forced by the coarse mesh boundary values, 

2) the interface boundary conditions, and 3) the phase dif- 
ferences in pressure oscillations. The u-component profile. 
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Fig. 15B, shows the reflection of energy from the east wall. 
Insufficient interface boundary conditions did not allow 
free exchange of mass and energy across the interface, 
thereby resulting in the reflections. 

In addition, Fig. 5 shows the pressure variations for 
the nested (one-way) fine mesh. These oscillations appear 
to follow the variations of the 1.5° uniform fine mesh 
(nested domain) spacing out to eight hours. After that time, 
the coarse mesh variations force the nested variations. This 
could lead to imbalances. 

When the two-way interaction scheme is used, the fine 
mesh has significant influence on the coarse mesh, Fig. 16. 

In one-way interaction the coarse mesh v-profile would show 
only a +.5 m/sec value for the three levels. It should be 
noted, though, that the variations of the fine mesh hori- 
zontal velocity components were not as large in the two-way 
as in the one-way. This was due to the transfer of energy 
out of the nested grid into the coarse grid. This allowed 
any oscillations present in the fine mesh region to propa- 
gate to the coarse mesh region. These oscillations, though, 
aggravated the problems at the coarse mesh boundary wall. 

From 12 hours and beyond, the oscillations in both grids 
become unrealistically large. 

The last series of experiments was computed using the 
nested, two-way interacting model with a 18° x 750 m mountain 
present. This allows three grid points in the coarse mesh 
and 11 grid points in the fine mesh to represent the mountain 
at each latitude. The resulting horizontal velocity profiles 
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after six hours for the fine and coarse meshes are shown in 
Figs. 17 and 18. The mountain has developed to half its 
height by this time. The v-profiles in Fig. 17A show the 
general upstream northward deflection and downstream south- 
ward deflection. In addition, the upstream flow appears to 
be deflected northward only in the immediate vicinity of the 
mountain. This is different than the coarse mesh solution 
which shows deflections three gird points upstream from the 
mountain. This narrowing of the velocity component response 
compares favorably with that of Hayes (1977) . One result 
which does not compare is the lack of dampening with height. 
This is probably a result of the improper v-component values 
from the coarse mesh being forced at the nested boundary in- 
terface . 

The u-profile (Fig. 17B) , for level 6(a = .92) does indi- 
cate an increase in the zonal flow over the mountain top, and 
then a return to the upstream value. The upper two levels 
(c = .42, .08), though, show an increase in the zonal flow 

downstream of the mountain. This is probably caused by the 
reflections from the east interface wall. 

The coarse mesh velocity profiles are shown in Figs. 18A 
and B. When compared to Figs. 9A and B, the one-way inter- 
acting mountain case, few significant differences are noticed. 
The v-profiles show the proper deflection north and southward, 
with the zonal flow increasing in the region of the mountain 
top. Since the flow in the region of the mountains is just 
an average of the nested region flow, the fine mesh solution 
must have been good for the first six hours. 
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In Figs. 19 and 20 are shown the v-component profiles 
and vertical velocities for the coarse mesh at 12 hours. 

The v-profiles show a large level l(a = .08) northward com- 
ponent with a compensating lower level 6(a = .92) southward 
component. The fine mesh v-profile is similar, except the 
deflections are spread wider over its domain. Part of the 
reason for the large upper level v-component is the inter- 
action of the mountain and the channel walls. As mentioned 
in section 1, large vertical velocities are created at the 
intersection with channel wall points, which causes an 
artificial meridionally oriented, vertical circulation. The 
same circulation apparently occurs in the fine mesh, due to 
interface, mountain intersection. This energy is then par- 
tially transferred to the coarse mesh via the two-way inter- 
action. The coarse mesh vertical velocities are shown in 
Fig. 20. The fine mesh influence is immediately apparent. 

By this time, the solutions to the nested grid model are 
dominated by unrealistic oscillations. Similar results were 
obtained in the one-way interaction solutions. 



40 



VI. SUMMARY AND CONCLUSIONS 



This thesis was undertaken to study the flow over moun- 
tains using a channeled model with an interior nested grid. 
The work involved two stages. The first was to develop an 
east-west channel model from a global model. The second was 
the nesting of a finer resolution grid. 

In channeling the modified Arakawa global model by 
Monaco and Williams (1975), oscillations and noise arose 
along the channel walls. The scheme C grid arrangement and 
flux form of the advective terms lead to certain boundary 
conditions. These proved to be adequate when no disturbance 
was present, but insufficient when either a mountain and/or 
nested grid were introduced. Difficulties also arose at the 
mountain- channel wall intersection. 

In nesting a finer resolution grid several things were 
learned. The numerical response in the ^-direction decreases 
with the finer resolution when compared to the coarse mesh 
solution. In addition, the two-way interaction does allow 
partial transfer of energy to the coarse mesh. Several prob- 
lems were also encountered. 

The grid arrangement and flux form of advective terms 
lead to similar boundary/interface conditions as in the 
coarse mesh. Difficulties did not arise at the interface 
walls due to extrapolation methods in providing boundary 
values, but they resulted due to not allowing free exchange 
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of energy and mass across the interface walls. This caused 
internal reflected fine mesh waves, which eventually de- 
stroyed the fine mesh solution. The mountain- interface 
wall intersection also accounted for some of the problems. 

As in the coarse mesh, large vertical velocities developed 
at these intersection points, which are partly responsible 
for the vertical meridional oscillations. 

Based on the above results, the following is recommended: 

(a) expand the channel in the north-south direction suffi- 
ciently to allow isolating the mountain. This would 
avoid the problems caused by the mountain channel wall 
intersection; 

(b) simplify the flux forms in the primitive equations for 
use over the whole domain; 

(c) modify the channel walls to consist of one latitude 
<j> , t , u row. This leads to an unequal predictive 
region for the v-component; 

(d) modify the interface walls in the nested grid to consist 
of one latitude (north and south walls) and one longi- 
tude (east and west walls) row. This leads to unequal 
predictive rows for v and u-components , respectively; 

(e) alter linear interpolation scheme to higher order scheme 
for providing the nested grid boundary with coarse mesh 
values . 

Although this thesis was not entirely successful in mathe- 
matically simulating flow over mountains, it did identify and 
detail some specific problems. It is hoped that this study 
will prove useful to those attempting a similar study in the 
future . 
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APPENDIX A 



NORTH AND SOUTH BOUNDARY CONSIDERATIONS 



The constraint of no net mass flux in the north-south 
direction is a primary factor to consider when determining 
boundary conditions for the north and south walls. Due to 
the cyclic condition in the east-west direction, the only 
equation of motion needed to be considered is equation (2-3) , 
the v-component equation: 



8 , TT v .3 / TTU \ , 3 ,TTV x , 3 ,1 TO x 

3t ( iK V) + 3C ( 1T V) + Tn^"m" V) + 3^ ( 5K v) 
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To assure no net mass flux, the above equation, integrated 
in £ and o , must balance. 

Following is an analysis of each term, and what must be 
done to accomplish the constraint condition. 

Term (1) : 



J f h<s v)dod5 ' It 



// 



@dodC - 0 



In the case where no mountain is present, starts at 

zero and should remain zero as long as the model initiali- 
zation is correct. When a mountain is present, there must 
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be a shift of mass to another region in order to satisfy 
the integral. 

Term (2) : 




fr^TT v)dod 5 - 0 



This is due to the cyclic condition east-west. 
Term (3) : 







/Also *** 1 



The first term on the right hand side is taken care of by 
reflecting the normal mass flux quantities along the bound- 
aries. In the no mountain case, the second term starts at 
zero and remains close to zero as the model advances in time. 
When a mountain is introduced, this term becomes non-zero. 
This problem is further accentuated by the imbalance in the 
pressure gradient term (7) . 

Term (4) : 

/ / io<= v)dod5 - ° 

5 a 

If d = 0 at a = 1 and a = 0 , this term is satisfied. 
Term (5): 

/ / 1ST dcdi + 0 

? a 
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This term is not equal to zero, but will primarily be bal- 
anced by the pressure gradient term. The term will further 
vary if a mountain is present. 

Term (6) : 
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The first term on the right hand side starts at zero (v=0) , 
and will remain close to zero as the model advances in time. 
The second term is not zero, because the integrated scaling 
factor is non-zero. The contribution, though, is very small 
compared to the other terms in the equation. Its importance, 
however, increases as the velocity components increase. 

Since the other terms are not exactly zero, this term adds 
to the imbalance. 

Term (7) : 

/ / E l H + 00 H ldodS + 0 

5 o 

This term is non-zero and primarily balanced by the £-compo- 
nent of velocity. In order to prevent mass variations, 
though, its initial balance along the boundaries must be 
maintained. The initial terrain pressure gradient was rela- 
tively simple to maintain but the initial geopotential height 
gradient was not so simple. 
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Term (3) : 





dad£ = 0 



This model is 



frictionless . 



APPENDIX B 



DIVERGENCE BOUNDARY CONDITION 

Free exchange of mass and energy across the nested grid 
interface is a necessity to prevent reflections from the 
walls destroying the fine mesh solutions. In addition to 
the two-way interaction scheme, a divergence boundary condi- 
tion after Madala (1977) was attempted. 

This boundary condition entailed calculating the coarse 
mesh divergence along the interface, and interpolating spa- 
tially and temporally to determine the fine mesh mass flux 
quantity on the outside of the interface. The following 
formulas were used for the east and south interfaces respec- 
tively : 



(iru — ), , = DIV + (ttu ^) t , , - (ttv ^) t ,,, 

v n I,J coarse v n'l-l.J ' m / I,J+l 
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The interface divergence boundary conditions were not applied 
along the north and west interfaces, due to the problems men- 
tioned in chapter IV. With this boundary condition used 
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along the east and south walls, results were still unsatis- 
factory. Reflections were still occurring, indicating a 
free transfer of mass and energy was not present. 

Madala (1977) apparently solved this problem by using 
weighted extrapolated fine mesh U and V values, along with 
coarse mesh values, in determining the mass flux divergence 
terms for the fine mesh boundary. This modified the forc- 
ing of the coarse mesh divergence fields at the interface. 
It should be noted that Madala used only the tt , T , and 
cp row of the scheme C grid for the boundary. In addition, 
he used a simplified flux form, thereby he did not need to 
specify as many external boundary points. This makes the 
north and west interface conditions exactly the same as the 
east and south, whereas in this research they would have 
been different. 
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Figure 1. The sigma (a) coordinate system as used 
in the 6 level model. P t is 200 mb. 
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Figure 2. 



Five possible grid distributions of the dependent 
variables $ , T , V , U (after Arakawa, 1974) . 
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PRESSURE 




Figure 3. Pressure variations of the coarse mesh. 

w/o (P-/YP) 

with (P-AP) 
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Figure 4. Temperature variations of the coarse mesh. 

w/o (P-AP) 

with (P-AP) 
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uniform fine mesh (nested domain) 
nested fine mesh 
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Figure 6B. South channel wall of coarse mesh. 
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Figure 7. Nested grid system. 



COARSE GRID 




Figure 8. 



Time step routines for a 45 minute sequence. 
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Figure 9A. V-components for coarse mesh with mountain A. 

time = 6 hours 




Figure 9B . U-components for coarse mesh with mountain A. 

time = 6 hours 

a = . 08 (solid line) a = ,42 (dashed line) 

a = .92 (dotted line) 
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Figure 10. Vertical velocity (j^r) profile for coarse mesh with 
mountain A at level 11 5 (o = .92). 
time = 6 hours 
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Figure 11. Temperature profiles for the six levels of 
the coarse mesh model with mountain A. 
time = 6 hours 
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Figure 12A. V-components for coarse mesh with mountain A. 

time = 12 hours 




Figure 12B. U-components for coarse mesh with mountain A. 

time = 12 hours 

a = .08 (solid line) a = .42 (dashed line) 
a = .92 (dotted line) 
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Figure 13. Vertical velocity (S profile for coarse mesh with 
mountain A at level 5 (a = .92). 
time = 12 hours 
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Figure 15A. 



- 10 *- 

V-components for fine mesh (one-way) with no 
mountain. 

time = 6 hours 




I 



Figure 15B. U-components for fine mesh (one-way) with no 
mountain. 

time = 6 hours 

a = .08 (solid line) a = .42 (dashed line) 

a = .92 (dotted line) 
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Figure 16. V- components for coarse mesh (two-way) with 
no mountain. 

time = 6 hours 

a = .08 (solid line) a = . 42 (dashed line) 

a = . 92 (dotted line) 
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Figure 17A. V- components for fine mesh (two-way) with 
mountain A. 

time = 6 hours 




Figure 17B. U-components for fine mesh (two-way) with 
mountain A. 

time = 6 hours 

a = .08 (solid line) a = .42 (dashed line) 

a = .92 (dotted line) 
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Figure 18A. 



V-components for coarse mesh (two-way) with 
mountain A. 

time = 6 hours 




Figure 18B. U-components for coarse mesh (two-way) with 
mountain A. 

time = 6 hours 

c = .08 (solid line) a = .42 (dashed line) 

a = .92 (dotted line) 
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Figure 19. V-components for coarse mesh (two-way) with 
mountain A. 

time = 12 hours 

a = . 08 (solid line) a = .42 (dashed line) 

o = .92 (dotted line) 
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Figure 20. Vertical velocity (p^r) profile for coarse mesh (two-way) 
with mountain A at level 5 (o = .92). 
time = 12 hours 
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